DESI DR1 Stream Tutorial¶

Welcome to the 2025 version of our stellar stream characterization tutorial notebook. This notebook will walk you through using data from the DESI Milky Way Survey data from the first full public release. For more information on the DESI DR1 stellar catalogue from the Milky Way Survey, see the README.md file.

Coments like the one below highlight the pieces of code you're expected to tweak in your search for stellar streams.

#-----------------------#

Example_variable = 1000000000000e10 #Look here for stuff to change throughout the notebook!

#-----------------------#

Below, we're going to import the packages we'll be using.

InĀ [1]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scipy as sp
import scipy.stats as stats
from astropy.io import fits
from astropy import table
import matplotlib
import matplotlib.patheffects as path_effects
import importlib
import stream_functions as stream_funcs
importlib.reload(stream_funcs)
import emcee
import corner
from astropy.table import Table
from astropy import units as u
from collections import OrderedDict
import time
from scipy import optimize, stats
import matplotlib.colors as mcolors
colors = mcolors.CSS4_COLORS
color_names = list(colors.keys())
import streamTutorial as st
import copy
importlib.reload(st);

Next, let's input our DESI MWS DR1 data.

InĀ [2]:
importlib.reload(st)
# Add the path to the DESI data and STREAMFINDER data below
#-----------------------#
desi_path = '../mwsall-pix-iron.fits'
sf_path = './data/streamfinder_gaiadr3.fits'
#-----------------------#

Data = st.Data(desi_path, sf_path)

print('Now our desi data has been loaded under Data.desi_data')


#      _______
#    (  ^   ^  )šŸ‘
#     |   ~   |
#     |_______|
#   ~ "Well done!"
Length of DESI Data before Cuts: 6372607
Length after NaN cut: 4075716
Now our desi data has been loaded under Data.desi_data

Pick Stream¶

We'll pick a stream to work with by entering it in

#-----------------------#
st.stream(Data, streamName='SoI-I21', streamNo=42)
#-----------------------#

Initializing this stream object will add an attribute to the Data object, accessed either through

Data.confirmed_sf_and_desi or SoI.data.confirmed_sf_and_desi

Here is the table to reference in picking your stream:

SF3 Stream Index Stream # SF3 Stars # SF3 Stars in DESI # of SF3 Stars with VHel Galstream Stream name Metallicities Distance_kpc
1 C-20 29.0 8.0 8.0 C-20-I24 -2.93
3 NGC 173.0 3.0 7.0 False -1.32
4 ATLAS 208.0 3.0 80.0 AAU-ATLAS -2.2
6 Kwando 138.0 32.0 62.0 Kwando-I21 -2.29 8.2
8 New-2 90.0 20.0 18.0 New-2-I24 -2.07
10 Gaia-12 46.0 9.0 5.0 Gaia-12-I21 -3.28 13.3
17 New-3 120.0 3.0 10.0 New-3-I24 -2.71
19 Leiptr 412.0 26.0 37.0 Leiptr-I21 -2.17 8.9
22 New-6 146.0 12.0 12.0 New-6-I24 -2.18
24 C-25 121.0 40.0 25.0 C-25-I24 -2.3
25 C-11 112.0 33.0 9.0 C-11-I24 -2.91
27 New-7 37.0 15.0 11.0 New-7-I24 -2.23
28 New-8 7.0 1.0 1.0 New-8-I24 -1.99
29 New-9 10.0 6.0 7.0 New-9-I24 -1.95
30 C-10 158.0 49.0 5.0 C-10-I24 -0.91
31 New-10 47.0 1.0 7.0 New-10-I24 -0.99
32 New-11 18.0 7.0 6.0 New-11-I24 -2.03
33 Gaia-10 141.0 45.0 26.0 Gaia-10-I21 -1.4 17.4
34 Gjoll 607.0 18.0 40.0 NGC3201-P21 -1.63
35 C-9 183.0 70.0 35.0 C-9-I24 -0.72
36 C-24 244.0 86.0 27.0 C-24-I24 -0.93
39 Slidr 330.0 148.0 34.0 Slidr-I21 -1.7 2.4
41 Ylgr 919.0 14.0 20.0 Ylgr-I21 -2.09 8.6
42 Sylgr 256.0 109.0 33.0 Sylgr-I21 -2.92 2.4
43 New-15 152.0 74.0 23.0 New-15-I24 -2.01
47 Fjorm 297.0 104.0 29.0 M68-P19 -2.23
48 Orphan 867.0 113.0 247.0 Orphan-K23 -1.9
49 Gaia-1 200.0 44.0 28.0 Gaia-1-I21 -1.8 5.0
50 C-23 29.0 12.0 6.0 C-23-I24 -2.36
51 LMS-1 358.0 89.0 29.0 LMS1 -2.09 16.2
52 C-22 39.0 17.0 11.0 C-22-I24 -2.05
53 GD-1 1468.0 670.0 323.0 GD-1-I21 -2.49 9.2
56 New-17 13.0 2.0 0.0 New-17-I24 -2.35
57 NGC5466 43.0 14.0 6.0 NGC5466-J21 -1.98 16.6
58 New-18 107.0 33.0 1.0 New-18-I24 -0.94
59 Gaia-6 145.0 63.0 14.0 Gaia-6-I21 -1.53 8.5
60 New-19 40.0 22.0 4.0 New-19-I24 -0.79
61 Pal-5 129.0 33.0 69.0 Pal5 -1.41
62 M5 91.0 33.0 11.0 M5-G19 -1.29 15.9
63 Kshir 141.0 50.0 7.0 Kshir-I21 -1.83 10.1
64 Svol 234.0 59.0 16.0 Svol-I21 -1.98 8.9
65 Gaia-9 233.0 78.0 25.0 Gaia-9-I21 -2.21 3.8
68 M92 202.0 140.0 23.0 M92-I21 -2.31 9.2
70 Gaia-11 84.0 11.0 3.0 Gaia-11-I21 -1.19 11.5
71 Hrid 666.0 46.0 29.0 Hrid-I21 -1.13 3.2
74 New-21 762.0 61.0 3.0 New-21-I24 -0.87
75 Phlegethon 632.0 85.0 41.0 Phlegethon-I21 -2.19 3.5
76 NGC7089 15.0 2.0 2.0 M2-I21 -1.65
79 New-23 44.0 4.0 6.0 New-23-I24 -2.25
80 New-24 56.0 1.0 7.0 New-24-I24 -0.64
82 New-26 43.0 11.0 4.0 New-26-I24 -0.17
84 NGC7492 26.0 1.0 1.0 NGC7492-I24 -1.78
85 C-19 46.0 6.0 12.0 C-19-I21 -3.58 18.0
86 Jhelum 986.0 1.0 160.0 Jhelum-a/b -2.12 10.9
InĀ [3]:
importlib.reload(stream_funcs)

#-----------------------#
streamName= 'Gaia-9-I21'#'Fjorm-I21' #'Sylgr-I21'
streamNo= 65# 47 #42
#-----------------------#

importlib.reload(st)
SoI = st.stream(Data, streamName, streamNo)
Importing galstreams module...
Initializing galstreams library from master_log... 
3.7774496996789693
Creating combined DataFrame of SF and DESI
No original data available for comparison - cut_confirmed_sf_and_desi is empty
Number of stars in SF: 233, Number of DESI and SF stars: 6
Saved merged DataFrame as self.data.confirmed_sf_and_desi
Stars only in SF3: 227

Let us take a look at our stellar stream below:

Figure 1a¶

InĀ [4]:
importlib.reload(st)
plt_soi = st.StreamPlotter(SoI)

plt_soi.on_sky(stream_frame=False)
/Users/nasserm/Documents/vscode/research/streamTut/DESI-DR1_streamTutorial/streamTutorial.py:762: UserWarning: You passed a edgecolor/edgecolors ('k') for an unfilled marker ('x').  Matplotlib is ignoring the edgecolor in favor of the facecolor.  This behavior may change in the future.
  ax.scatter(

No description has been provided for this image

We also calculated the stream coordinates $\phi_1$ and $\phi_2$. This is achieved by rotating the right ascension and declination such that the length of the stream lies along $\phi_2 \sim 0$, and the center of the stream is $\phi_1 \sim 0$.

Figure 1b¶

InĀ [5]:
plt_soi.on_sky(stream_frame=True)
No description has been provided for this image

Great! Lets move on now.

We have lots of DESI data, lets cut that down with the following steps:

  • Perform an RA/DEC cut to remove stars on different parts of the sky
  • Perform a distance cut to remove stars too nearby
  • Trim stars that are too metal rich
  • Trim stars that fall outside of the best-fit stellar population model
  • Trim stars with kinematics nothing like the STREAMFINDER stars

Spatial Cuts¶

The first two bullet points can be achieved using a handy function we have stored in stream_functions.py

stream_funcs.threeD_max_min_mask.

InĀ [6]:
importlib.reload(st)
selection_fine = st.Selection(SoI.data.desi_data)

# We want to get the ra and dec from STREAMFINDER stars for this function so we can cut around it

#-----------------------#
ra_cut = 5 #deg
dec_cut = 5 #deg
#-----------------------#

selection_fine.add_mask(name='3D',
        mask_func=lambda df: stream_funcs.threeD_max_min_mask(
        df['TARGET_RA'],        
        df['TARGET_DEC'],         
        df['PARALLAX'],       
        df['PARALLAX_ERROR'],   
        SoI.data.SoI_streamfinder['RAdeg'],         
        SoI.data.SoI_streamfinder['DEdeg'],        
        SoI.min_dist,              
        ra_cut,dec_cut) #<-----------------------# dec cut, ra cut based on streamfinder star values
)
Selection object created for DataFrame with 4075716 rows.
Mask added: '3D'

Once we've done all our cuts, we can get all our masks in one using mask = selection_fine.get_masks(['mask 1', ... 'mask n'])

Figure 2¶

InĀ [7]:
importlib.reload(st)
# Lets take a look at our trimmed DESI data

mask = selection_fine.get_masks(['3D'])

# Using the new method that replaces the 4-line pattern
trimmed_stream = SoI.mask_stream(mask)

plt_trim = st.StreamPlotter(trimmed_stream)
plt_trim.on_sky(showStream=True, background=True, stream_frame=False)
plt_trim.on_sky(showStream=True, background=True)
...'3D' selected 105268 stars
Selection for specified masks: 105268 / 4075716 stars.
No stars were cut - cut_confirmed_sf_and_desi is empty
Number of stars in SF: 233, Number of DESI and SF stars: 6
Saved merged DataFrame as self.data.confirmed_sf_and_desi_b
Stars only in SF3: 227
/Users/nasserm/Documents/vscode/research/streamTut/DESI-DR1_streamTutorial/streamTutorial.py:762: UserWarning: You passed a edgecolor/edgecolors ('k') for an unfilled marker ('x').  Matplotlib is ignoring the edgecolor in favor of the facecolor.  This behavior may change in the future.
  ax.scatter(

No description has been provided for this image
No description has been provided for this image

Lets also look at our parallax cut. We use the calculation of plx - 2 * plx_err < 1/min_dist. This simply gets rid of nearby stars that, after subtracting 2 times the parallax error, are still below minimum parallax.

Figure 3¶

InĀ [8]:
plt_trim.plx_cut()
/Users/nasserm/Documents/vscode/research/streamTut/DESI-DR1_streamTutorial/streamTutorial.py:828: UserWarning: You passed a edgecolor/edgecolors ('k') for an unfilled marker ('x').  Matplotlib is ignoring the edgecolor in favor of the facecolor.  This behavior may change in the future.
  ax.scatter(

No description has been provided for this image

Now that we've visualized the stream in the stream-frame, we can get a bit more specific with our on-sky cuts. We can narrow in by only keeping stars within some range of $\phi_2$ = 0.

InĀ [13]:
#-----------------------#
phi2_wiggle = 10 #[deg]
#-----------------------#

selection_fine.add_mask(name='phi2',
        mask_func=lambda df: (df['phi2'] < phi2_wiggle) & (df['phi2'] > -1*phi2_wiggle))
Mask added: 'phi2'

Figure 4¶

InĀ [14]:
importlib.reload(st)
# Lets take a look at our trimmed DESI data
mask = selection_fine.get_masks(['3D', 'phi2'])

trimmed_stream = SoI.mask_stream(mask)

plt_trim = st.StreamPlotter(trimmed_stream)
plt_trim.on_sky(showStream=True, background=True, stream_frame=False)
plt_trim.on_sky(showStream=True, background=True)
...'3D' selected 105268 stars
...'phi2' selected 721628 stars
Selection for specified masks: 38899 / 4075716 stars.
No stars were cut - cut_confirmed_sf_and_desi is empty
Number of stars in SF: 233, Number of DESI and SF stars: 6
Saved merged DataFrame as self.data.confirmed_sf_and_desi_b
Stars only in SF3: 227
/Users/nasserm/Documents/vscode/research/streamTut/DESI-DR1_streamTutorial/streamTutorial.py:762: UserWarning: You passed a edgecolor/edgecolors ('k') for an unfilled marker ('x').  Matplotlib is ignoring the edgecolor in favor of the facecolor.  This behavior may change in the future.
  ax.scatter(

No description has been provided for this image
No description has been provided for this image

Metallicity Cut¶

Stars in a stellar stream all share some chemical characteristics, and we can use this to cut down on the contaminents in our field. Lets grab the stream's known metallicity from STREAMFINDER and make a cut using that information.

InĀ [15]:
sf3_table = pd.read_csv('./data/sf3_only_table.csv')


sf_streamname = streamName.rsplit('-', 1)[0] #If this fails, manually enter the stream name from the first table.
metallicity = sf3_table[sf3_table['Stream'] == sf_streamname]['Metallicities'].values
print(f"{sf_streamname} Metallicity: {metallicity}")
print(f'Mass Fraction (Z) Guess: {0.0181 * 10 ** metallicity}')
Gaia-9 Metallicity: [-2.21]
Mass Fraction (Z) Guess: [0.0001116]

Figure 5¶

InĀ [16]:
plt_trim.plot_params['background']['alpha']=0.05
plot = plt_trim.feh_plot(showStream=True, background=True)
ax = plot[1]
ax.axhline(metallicity, color='g', linestyle='solid', label='Metallicity from SF3', alpha=0.5)
ax.legend()
/Users/nasserm/Documents/vscode/research/streamTut/DESI-DR1_streamTutorial/streamTutorial.py:1006: UserWarning: You passed a edgecolor/edgecolors ('k') for an unfilled marker ('x').  Matplotlib is ignoring the edgecolor in favor of the facecolor.  This behavior may change in the future.
  ax.scatter(

Out[16]:
<matplotlib.legend.Legend at 0x17781d890>
No description has been provided for this image

Now lets trim stars that are too metal rich!

InĀ [45]:
#-----------------------#
feh_cut = -1 # [dex]
#-----------------------#

selection_fine.add_mask(name='feh',
        mask_func=lambda df: (df['FEH'] < feh_cut))
Mask added: 'feh'

Figure 6¶

InĀ [46]:
#del trimmed_stream, plt_trim # clear memory
importlib.reload(st)
mask = selection_fine.get_masks(['3D', 'phi2', 'feh'])

trimmed_stream = SoI.mask_stream(mask)


plt_trim = st.StreamPlotter(trimmed_stream)
plot = plt_trim.feh_plot(showStream=True, background=True)
ax = plot[1]
ax.axhline(metallicity, color='g', linestyle='solid', label='Metallicity from SF3', alpha=0.5)
ax.axhline(feh_cut, color='red', linestyle='--', label='Cut')
ax.legend(loc='upper right')
...'3D' selected 105268 stars
...'phi2' selected 721628 stars
...'feh' selected 866850 stars
Selection for specified masks: 16099 / 4075716 stars.
No stars were cut - cut_confirmed_sf_and_desi is empty
Number of stars in SF: 233, Number of DESI and SF stars: 6
Saved merged DataFrame as self.data.confirmed_sf_and_desi_b
Stars only in SF3: 227
/Users/nasserm/Documents/vscode/research/streamTut/DESI-DR1_streamTutorial/streamTutorial.py:1006: UserWarning: You passed a edgecolor/edgecolors ('k') for an unfilled marker ('x').  Matplotlib is ignoring the edgecolor in favor of the facecolor.  This behavior may change in the future.
  ax.scatter(

Out[46]:
<matplotlib.legend.Legend at 0x17e1f9f10>
No description has been provided for this image

Isochrone Cut¶

When searching for stellar streams, an isochrone cut helps isolate stars that are likely part of the same stellar population. Streams originate from disrupted star clusters or dwarf galaxies, whose stars share a common age and chemical composition. By selecting only stars lying close to a theoretical stellar evolution isochrone (in colour–magnitude space), we suppress contamination from unrelated field stars with different ages and metallicities.

The position of an isochrone in colour-magnitude space depends on the age of the stellar population (~shifts it up and down) and the metallicity (~shifts it left and right)

InĀ [51]:
#-----------------------#
colour_wiggle = 0.2 # [dex]
age = 10.5 # [Gyr] -> of form 10.0, 10.1, 10.2, etc

metallicity = -1.5 #metallicity # may want to adjust the metallicity of the isochrone manually
#-----------------------#
SoI.isochrone(metallicity, age=age) 

selection_fine.add_mask(name='iso',
        mask_func=lambda df: (stream_funcs.betw(SoI.data.desi_colour_idx, SoI.isochrone_fit(SoI.data.desi_abs_mag), colour_wiggle)))
Mass Fraction (Z): 0.0005723722564904767

using ./data/dotter/iso_a10.5_z0.00057.dat
Using distance gradient
/Users/nasserm/Documents/vscode/research/streamTut/DESI-DR1_streamTutorial/stream_functions.py:635: RuntimeWarning: invalid value encountered in log10
  R_mag = r_mag - 5 * np.log10(distances) + 5

Mask added: 'iso'

Figure 7¶

InĀ [52]:
importlib.reload(st)
mask = selection_fine.get_masks(['3D', 'phi2', 'feh', 'iso'])

trimmed_stream = SoI.mask_stream(mask)
trimmed_stream.isochrone(metallicity, age=age) 

plt_trim = st.StreamPlotter(trimmed_stream)
plt_trim.plot_params['background']['alpha']=0.1
plt_trim.iso_plot(wiggle=colour_wiggle, background=True, showStream=True)
...'3D' selected 105268 stars
...'phi2' selected 721628 stars
...'feh' selected 866850 stars
...'iso' selected 1047795 stars
Selection for specified masks: 12170 / 4075716 stars.
No stars were cut - cut_confirmed_sf_and_desi is empty
Number of stars in SF: 233, Number of DESI and SF stars: 6
Saved merged DataFrame as self.data.confirmed_sf_and_desi_b
Stars only in SF3: 227
Mass Fraction (Z): 0.0005723722564904767

using ./data/dotter/iso_a10.5_z0.00057.dat
Using distance gradient
/Users/nasserm/Documents/vscode/research/streamTut/DESI-DR1_streamTutorial/streamTutorial.py:1060: UserWarning: You passed a edgecolor/edgecolors ('k') for an unfilled marker ('x').  Matplotlib is ignoring the edgecolor in favor of the facecolor.  This behavior may change in the future.
  ax.scatter(

No description has been provided for this image

Kinematic Cuts¶

DESI's radial velocity measurements are generally more reliable than from the Gaia RVS data (which the STREAMFINDER Catalogue is based off of). We may see some obvious outliers in the velocity data, and we will want to cut those out.

Figure 8¶

InĀ [53]:
importlib.reload(st)
plt_trim.plot_params['background']['alpha']=0.08
plt_trim.kin_plot(showStream=True, background=True, show_sf_only=False)
/Users/nasserm/Documents/vscode/research/streamTut/DESI-DR1_streamTutorial/streamTutorial.py:885: UserWarning: You passed a edgecolor/edgecolors ('k') for an unfilled marker ('x').  Matplotlib is ignoring the edgecolor in favor of the facecolor.  This behavior may change in the future.
  ax[0].scatter(

/Users/nasserm/Documents/vscode/research/streamTut/DESI-DR1_streamTutorial/streamTutorial.py:890: UserWarning: You passed a edgecolor/edgecolors ('k') for an unfilled marker ('x').  Matplotlib is ignoring the edgecolor in favor of the facecolor.  This behavior may change in the future.
  ax[1].scatter(

/Users/nasserm/Documents/vscode/research/streamTut/DESI-DR1_streamTutorial/streamTutorial.py:895: UserWarning: You passed a edgecolor/edgecolors ('k') for an unfilled marker ('x').  Matplotlib is ignoring the edgecolor in favor of the facecolor.  This behavior may change in the future.
  ax[2].scatter(

No description has been provided for this image
InĀ [60]:
#-----------------------#
vgsr_max=-100 #[km/s]
vgsr_min=-300
#-----------------------#

selection_fine.add_mask(name='VGSR',
        mask_func=lambda df: (df['VGSR'] < vgsr_max) & (df['VGSR'] > vgsr_min))
Mask added: 'VGSR'
InĀ [61]:
#-----------------------#
pmra_wiggle = 6 # [mas/yr]
#-----------------------#

selection_fine.add_mask(
    name='PMRA',
    mask_func=lambda df: (
        (
            df['PMRA'] >= np.interp(df['phi1'], [SoI.data.SoI_streamfinder['phi1'].min(), SoI.data.SoI_streamfinder['phi1'].max()],
                                            [SoI.data.SoI_streamfinder.loc[SoI.data.SoI_streamfinder['phi1'].idxmin(), 'pmRA'] - pmra_wiggle,
                                             SoI.data.SoI_streamfinder.loc[SoI.data.SoI_streamfinder['phi1'].idxmax(), 'pmRA'] - pmra_wiggle])
        ) &
        (
            df['PMRA'] <= np.interp(df['phi1'], [SoI.data.SoI_streamfinder['phi1'].min(), SoI.data.SoI_streamfinder['phi1'].max()],
                                            [SoI.data.SoI_streamfinder.loc[SoI.data.SoI_streamfinder['phi1'].idxmin(), 'pmRA'] + pmra_wiggle,
                                             SoI.data.SoI_streamfinder.loc[SoI.data.SoI_streamfinder['phi1'].idxmax(), 'pmRA'] + pmra_wiggle])
        ) ))
Mask added: 'PMRA'
InĀ [62]:
#-----------------------#
pmdec_wiggle = 6 # [mas/yr]
#-----------------------#


selection_fine.add_mask(
    name='PMDEC',
    mask_func=lambda df: (
        (
            df['PMDEC'] >= np.interp(df['phi1'], [SoI.data.SoI_streamfinder['phi1'].min(), SoI.data.SoI_streamfinder['phi1'].max()],
                                            [SoI.data.SoI_streamfinder.loc[SoI.data.SoI_streamfinder['phi1'].idxmin(), 'pmDE'] - pmdec_wiggle,
                                             SoI.data.SoI_streamfinder.loc[SoI.data.SoI_streamfinder['phi1'].idxmax(), 'pmDE'] - pmdec_wiggle])
        ) &
        (
            df['PMDEC'] <= np.interp(df['phi1'], [SoI.data.SoI_streamfinder['phi1'].min(), SoI.data.SoI_streamfinder['phi1'].max()],
                                            [SoI.data.SoI_streamfinder.loc[SoI.data.SoI_streamfinder['phi1'].idxmin(), 'pmDE'] + pmdec_wiggle,
                                             SoI.data.SoI_streamfinder.loc[SoI.data.SoI_streamfinder['phi1'].idxmax(), 'pmDE'] + pmdec_wiggle])
        ) ))
Mask added: 'PMDEC'

It is sometimes easier to do wider cuts in PMRA/PMDEC for finding the stream by-eye, you can use those masks instead below.

InĀ [63]:
pmra_max = -2 # [mas/yr]
pmra_min = -12

pmdec_max = 10 # [mas/yr]
pmdec_min = 0


#-----------------------#

selection_fine.add_mask(name='PMRA_wide',
        mask_func=lambda df: (df['PMRA'] < pmra_max) & (df['PMRA'] > pmra_min))

selection_fine.add_mask(name='PMDEC_wide',
        mask_func=lambda df: (df['PMDEC'] < pmdec_max) & (df['PMDEC'] > pmdec_min))
Mask added: 'PMRA_wide'
Mask added: 'PMDEC_wide'

Lets look at our cuts!¶

InĀ [64]:
selection_fine.list_masks()
Active masks:
- 3D
- phi2
- feh
- iso
- VGSR
- PMDEC
- PMRA_wide
- PMDEC_wide
- PMRA

Figure 9¶

InĀ [67]:
importlib.reload(st)
# Lets take a look at our trimmed DESI data
mask = selection_fine.get_masks(['3D', 'phi2', 'feh','iso', 'VGSR', 'PMRA_wide', 'PMDEC_wide']) # TIP 'PMRA_wide', 'PMDEC_wide' can be used instead if defined

# Using the new method that replaces the 4-line pattern
# VGSR is now automatically computed for confirmed_sf_not_desi within mask_stream()
trimmed_stream = SoI.mask_stream(mask)

plt_trim = st.StreamPlotter(trimmed_stream)
plt_trim.plot_params['sf_in_desi']['alpha']= 0.9
plt_trim.plot_params['background']['alpha']= 0.9
plt_trim.plot_params['background']['s'] = 10
#-----------------------#
plt_trim.kin_plot(showStream=False, show_sf_only=False, background=True) # You can change showStream to True to see where the STREAMFINDER stars are in this space.
                                                                       # If there aren't many SF stars in DESI, you can show the stars not in DESI as a guide by setting show_sf_only to True
#-----------------------#
...'3D' selected 105268 stars
...'phi2' selected 721628 stars
...'feh' selected 866850 stars
...'iso' selected 1047795 stars
...'VGSR' selected 653284 stars
...'PMRA_wide' selected 1625961 stars
...'PMDEC_wide' selected 715060 stars
Selection for specified masks: 69 / 4075716 stars.
No stars were cut - cut_confirmed_sf_and_desi is empty
Number of stars in SF: 233, Number of DESI and SF stars: 6
Saved merged DataFrame as self.data.confirmed_sf_and_desi_b
Stars only in SF3: 227
No description has been provided for this image

Lets look at an example of a stream (Sylgr-I21) being spotted once our initial cuts have been made. We can see lines of overdensities in all three kinematic dimensions!

Furthermore, we can overlay the STREAMFINDER stars to see that this overdensity lies on the previously disovered stream, and which SF stars we have cutout.

Stream Image Stream Image

If you don't see a stream, there is no point in continuing

These cuts are quite restrictive for the sake of spotting a stream by-eye. We'll widen them now to have more stars to run our mixture model on.

InĀ [68]:
import copy
selection_mcmc = copy.copy(selection_fine)
InĀ [69]:
## Make your cuts wider than the ones above.
#-----------------------#
colour_wiggle = 0.2 # [dex]

vgsr_max=  -100 # [km/s]
vgsr_min= -400

pmra_max = 0 # [mas/yr]
pmra_min = -30

pmdec_max = 20 # [mas/yr]
pmdec_min = -35

feh_wide = -1
#-----------------------#

selection_fine.add_mask(name='feh_wide',
        mask_func=lambda df: (df['FEH'] < feh_wide))

selection_mcmc.add_mask(name='iso_wide',
        mask_func=lambda df: (stream_funcs.betw(SoI.data.desi_colour_idx, SoI.isochrone_fit(SoI.data.desi_abs_mag), colour_wiggle)))

selection_mcmc.add_mask(name='VGSR_wide',
        mask_func=lambda df: (df['VGSR'] < vgsr_max) & (df['VGSR'] > vgsr_min))

selection_mcmc.add_mask(name='PMRA_wide',
        mask_func=lambda df: (df['PMRA'] < pmra_max) & (df['PMRA'] > pmra_min))

selection_mcmc.add_mask(name='PMDEC_wide',
        mask_func=lambda df: (df['PMDEC'] < pmdec_max) & (df['PMDEC'] > pmdec_min))
Mask added: 'feh_wide'
Mask added: 'iso_wide'
Mask added: 'VGSR_wide'
Mask added: 'PMRA_wide'
Mask added: 'PMDEC_wide'
InĀ [70]:
mcmc_mask = selection_mcmc.get_masks(['3D', 'phi2', 'feh_wide','iso_wide', 'VGSR_wide', 'PMRA_wide', 'PMDEC_wide'])
...'3D' selected 105268 stars
...'phi2' selected 721628 stars
...'feh_wide' selected 866850 stars
...'iso_wide' selected 1047795 stars
...'VGSR_wide' selected 655452 stars
...'PMRA_wide' selected 2685256 stars
...'PMDEC_wide' selected 3989672 stars
Selection for specified masks: 1666 / 4075716 stars.

Below, lets look at our final selection of stars before we run our optimization.

Figure 10¶

InĀ [71]:
importlib.reload(st)
tomcmc_stream = SoI.mask_stream(mcmc_mask)

plt_trim = st.StreamPlotter(tomcmc_stream)
plt_trim.plot_params['sf_in_desi']['alpha']= 0.9
plt_trim.plot_params['background']['alpha']= 0.4
plt_trim.plot_params['background']['s'] = 1
#-----------------------#
plt_trim.sixD_plot(showStream=True, show_sf_only=False, background=True) # You can change showStream to True to see where the STREAMFINDER stars are in this space.
                                                                       # If there aren't many SF stars in DESI, you can show the stars not in DESI as a guide by setting show_sf_only to True
#-----------------------#
No stars were cut - cut_confirmed_sf_and_desi is empty
Number of stars in SF: 233, Number of DESI and SF stars: 6
Saved merged DataFrame as self.data.confirmed_sf_and_desi_b
Stars only in SF3: 227
Out[71]:
(<Figure size 1000x1500 with 5 Axes>,
 array([<Axes: ylabel='$\\phi_2$'>, <Axes: ylabel='V$_{GSR}$ (km/s)'>,
        <Axes: ylabel='$\\mu_{\\alpha}$ [mas/yr]'>,
        <Axes: ylabel='$\\mu_{\\delta}$ [mas/yr]'>,
        <Axes: xlabel='$\\phi_1$', ylabel='[Fe/H]'>], dtype=object))
No description has been provided for this image

Optimizing¶

If you see a stream above, you may continue with the rest of the notebook!

Before jumping into running MCMC, we will find a good starting point using scipy.minimize. This reduces the amount of time we'll need to run the MCMC, and generally means we'll get a better fit.

Streams aren't flat lines in $\phi_1$ dynamically, so we want to allow our fits to vary along the track.

#-----------------------#
no_of_spline_points = 5
#-----------------------#

Typically we want between 3 and 5 spline points.

InĀ [100]:
importlib.reload(st)
#-----------------------#
no_of_spline_points = 5

# Optional: Set custom spline range (if not specified, uses StreamFinder data range)
phi1_min = 2  # Set to a specific value to override the leftmost spline location
phi1_max = 28  # Set to a specific value to override the rightmost spline location

#-----------------------#

# Define truncation parameters based on our selection cuts
truncation_params = {
    'vgsr_min': vgsr_min, 'vgsr_max': vgsr_max,
    'feh_min': -4.0, 'feh_max': feh_wide,
    'pmra_min': pmra_min, 'pmra_max': pmra_max,
    'pmdec_min': pmdec_min, 'pmdec_max': pmdec_max
}

MCMeta = st.MCMeta(no_of_spline_points, tomcmc_stream, trimmed_stream.data.confirmed_sf_and_desi, 
                   truncation_params=truncation_params, phi1_min=phi1_min, phi1_max=phi1_max)
Making stream initial guess based on galstream and STREAMFINDER...
Stream VGSR dispersion from trimmed SF: 5.45 km/s
Stream mean metallicity from trimmed SF: -2.10 +- 0.196 dex
Stream PMRA dispersion from trimmed SF: 0.88 mas/yr
Stream PMDEC dispersion from trimmed SF: 1.33 mas/yr
Making background initial guess...
Background velocity: -155.39 +- 49.73 km/s
Background metallicity: -1.71 +- 0.456 dex
Background PMRA: -4.36 +- 3.89 mas/yr
Background PMDEC: -5.72 +- 4.54 mas/yr

We're doing a mixture model, but what are we mixing? Lets take a look at our intial guesses below:

Figure 11¶

InĀ [101]:
plt_mcmeta = st.StreamPlotter(MCMeta)

# Show our initial parameter guesses
print("Initial Parameters:")
for param, value in MCMeta.initial_params.items():
    if isinstance(value, np.ndarray):
        print(f"{param}: {value}")
    else:
        print(f"{param}: {value:.4f}")

# Plot the Gaussian mixture model based on our initial guesses
plt_mcmeta.gaussian_mixture_plot(showStream=True, background=True)
Initial Parameters:
lsigvgsr: 0.7363
vgsr_spline_points: [-227.97172772 -217.91524139 -221.40285946 -238.43458195 -269.01040885]
feh1: -2.1028
lsigfeh: -0.7085
lsigpmra: -0.0531
pmra_spline_points: [-7.14078642 -9.96807727 -9.64118496 -6.16010948  0.47514916]
lsigpmdec: 0.1239
pmdec_spline_points: [9.40057652 5.65807403 3.45369257 2.78743211 3.65929266]
bv: -155.3923
lsigbv: 1.6966
bfeh: -1.7062
lsigbfeh: -0.3410
bpmra: -4.3650
lsigbpmra: 0.5897
bpmdec: -5.7244
lsigbpmdec: 0.6572
Out[101]:
(<Figure size 900x900 with 4 Axes>,
 array([[<Axes: xlabel='V$_{GSR}$ (km/s)'>, <Axes: xlabel='[Fe/H]'>],
        [<Axes: xlabel='$\\mu_{RA}$ (mas/yr)'>,
         <Axes: xlabel='$\\mu_{DEC}$ (mas/yr)'>]], dtype=object))
No description has been provided for this image

You can see our initial guess above. showing the mixture of the stream model (blue) and the background model (orange) combining to make the overall distribution (black). Depending on your cuts, you may also notice that the gaussians are truncated for some parameters—this is done to allow us to cut down on the background stars while still modelling each population with some gaussian mean and spread.

Priors¶

An important step for our analysis is to provide strong priors about the background and stream. We will base those around our initial guess, using the STREAMFINDER catalogue to inform our priors.

InĀ [102]:
importlib.reload(st)
vgsr_range_wiggle = 10
pmra_range_wiggle = 5
pmdec_range_wiggle = 10
lsigvgsr_r = (-4, 4)

vgsr_ranges = [(v - vgsr_range_wiggle, v + vgsr_range_wiggle) for v in MCMeta.initial_params['vgsr_spline_points']]
pmra_ranges = [(v - pmra_range_wiggle, v + pmra_range_wiggle) for v in MCMeta.initial_params['pmra_spline_points']]
pmdec_ranges = [(v - pmdec_range_wiggle, v + pmdec_range_wiggle) for v in MCMeta.initial_params['pmdec_spline_points']]

#-----------------------#
feh1_r = (-3,-2)
#-----------------------#
lsigfeh_r = (-2,4)
lsigpmra_r = (-5, 1)
lsigpmdec_r = (-5, 1)
bv_r = (-400, 400)
lsigbv_r = (-2,4)
bfeh_r = (-4, 4)
lsigbfeh_r = (-2,2)
bpmra_r = (-20, 20)
lsigbpmra_r = (-2,3)
bpmdec_r = (-20, 10)
lsigbpmdec_r = (-2,3)

pstream_r = (0.005, 1.0)

prior = [
    pstream_r,              # pstream (single value)
    *vgsr_ranges,           # vgsr_spline_points (multiple values)
    lsigvgsr_r,             # lsigvgsr (single value)
    feh1_r, lsigfeh_r,      # feh parameters
    *pmra_ranges,           # pmra_spline_points (multiple values)  
    lsigpmra_r,             # lsigpmra (single value)
    *pmdec_ranges,          # pmdec_spline_points (multiple values)
    lsigpmdec_r,            # lsigpmdec (single value)
    bv_r, lsigbv_r, bfeh_r, lsigbfeh_r,  # background parameters
    bpmra_r, lsigbpmra_r, bpmdec_r, lsigbpmdec_r
]

MCMeta.priors(prior)

Now that we have our priors set, lets try running a basic minimization to find a good starting point for the MCMC. This will run from 10 to 60 seconds...

InĀ [103]:
MCMeta.scipy_optimize()
Running optimization...
Maximum number of function evaluations has been exceeded.

Optimized Parameters:
pstream: 0.0115
vgsr_spline_points: [-223.81793425 -221.61821854 -230.43612139 -228.75638726 -272.31749216]
sigvgsr: 24.5321
feh1: -2.0070
sigfeh: 0.2826
pmra_spline_points: [ -2.18778749 -11.19179392  -8.56063685  -7.17247656   0.60860223]
sigpmra: 0.9195
pmdec_spline_points: [19.31499566  4.20760611  4.04286177  6.85089452  6.47009953]
sigpmdec: 2.5767
bv: 175.3628
sigbv: 143.9422
bfeh: -1.4949
sigbfeh: 0.5687
bpmra: 2.8947
sigbpmra: 6.8304
bpmdec: -5.8216
sigbpmdec: 4.4956

Figure 13¶

InĀ [104]:
# First reload the module to get the fixed version
importlib.reload(st)


plt_optimized = st.StreamPlotter(MCMeta)
original_params = MCMeta.initial_params.copy()
optimized_for_plotting = MCMeta.sp_output.copy()

# Update MCMeta with optimized parameters for plotting
MCMeta.initial_params = optimized_for_plotting

print("\nPlotting optimized Gaussian mixture...")
# This should now work without the broadcasting error
plt_optimized.gaussian_mixture_plot(showStream=False, background=True)

# Restore original initial_params
MCMeta.initial_params = original_params

plt.suptitle('Optimized Gaussian Mixture Model', fontsize=14, y=0.98)
plt.tight_layout()
plt.show()

MCMeta.prior_validation()
Plotting optimized Gaussian mixture...
No description has been provided for this image
Your prior is good, you've found something!
All 70 walkers initialized successfully!

The above cell automatically checks if your prior was any good. If you get a message like the following:

prior_error.png

Try going back and adjusting the mentioned parameter.

Figure 14¶

InĀ [105]:
importlib.reload(st)

plt_enhanced = st.StreamPlotter(MCMeta) 
plt_enhanced.plot_params['sf_in_desi']['alpha'] = 0.9
plt_enhanced.plot_params['background']['alpha'] = 0.4
plt_enhanced.plot_params['background']['s'] = 1

fig, ax = plt_enhanced.sixD_plot(
    showStream=True, 
    show_sf_only=False, 
    background=True,
    show_initial_splines=True,
    show_optimized_splines=True,
    show_sf_errors=True  # Now enable error bars
)

plt.tight_layout()
plt.show()
No description has been provided for this image

EMCEE¶

Once we've optimized our intial guess with Scipy, we can run our MCMC sampler.

Only run this is you are interested in the MCMC sampler. This will take up to 20 minutes unless you reduce the number of burnin's or steps taken.

InĀ [106]:
importlib.reload(st)
from datetime import datetime
import os
note = 'test'
current_date = datetime.now().strftime("%y%m%d")
output_dir = f"./runs/{streamName}_{current_date}_{note}"
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

mcmc = st.MCMC(MCMeta, output_dir)
InĀ [107]:
%%time 
#-----------------------#
nproc = 32
nburnin = 5000
nstep = 5000
use_optimized_start = True  # Toggle: True for optimized results, False for initial guess
#-----------------------#

mcmc.run(nproc=nproc, nburnin=nburnin, nstep=nstep, use_optimized_start=use_optimized_start)
Using optimized parameters as starting positions...
Running burn-in with 5000 iterations. starting from optimized parameters...
Clipping all walker positions to be within prior ranges...
Took 95.7 seconds (1.6 minutes)
Now sampling with 5000 iterations
Took 87.0 seconds (1.5 minutes)
Getting flatchain...
CPU times: user 1min 53s, sys: 58.2 s, total: 2min 51s
Wall time: 3min 3s

Figure 15¶

InĀ [108]:
mcmc.show_chains()
No description has been provided for this image

Figure 16¶

InĀ [109]:
mcmc.show_corner()
WARNING:root:Too few points to create valid contours
No description has been provided for this image
InĀ [110]:
mcmc.print_result()
29
OrderedDict([('pstream', 0.012145763667586056), ('vgsr1', -226.91896470928137), ('vgsr2', -217.2252196217937), ('vgsr3', -227.41981876494893), ('vgsr4', -233.8826551822452), ('vgsr5', -263.047122948714), ('lsigvgsr', 0.9249192044680123), ('feh1', -2.0433849723565793), ('lsigfeh', -1.3518705705625096), ('pmra1', -6.958775801114271), ('pmra2', -9.978454030873168), ('pmra3', -9.015263144910275), ('pmra4', -7.298067695843081), ('pmra5', -4.155569655248742), ('lsigpmra', -0.0826696880854787), ('pmdec1', 9.2855245980585), ('pmdec2', 5.31564621924706), ('pmdec3', 3.105916786774497), ('pmdec4', 2.512317435240477), ('pmdec5', 1.1327716314762475), ('lsigpmdec', -0.1444537349135165), ('bv', 209.12913719714743), ('lsigbv', 2.176917709413055), ('bfeh', -1.4845591731953998), ('lsigbfeh', -0.24404068929637918), ('bpmra', 16.592038033174013), ('lsigbpmra', 1.0112094442352282), ('bpmdec', -5.829635415611426), ('lsigbpmdec', 0.6502205681528926)])
param             med         em         ep   exp(med)    exp(em)    exp(ep)
--------------------------------------------------------------------------------------
pstream         0.012     -0.003      0.003
vgsr1        -226.919     -7.796      6.506
vgsr2        -217.225     -7.068      5.629
vgsr3        -227.420     -2.679      3.377
vgsr4        -233.883     -3.265      2.860
vgsr5        -263.047     -7.240      3.039
lsigvgsr        0.925     -0.158      0.120      8.412     -2.569      2.678
feh1           -2.043     -0.030      0.025
lsigfeh        -1.352     -0.448      0.397      0.044     -0.029      0.066
pmra1          -6.959     -3.274      3.480
pmra2          -9.978     -1.090      0.988
pmra3          -9.015     -0.507      0.562
pmra4          -7.298     -0.336      0.348
pmra5          -4.156     -0.274      0.811
lsigpmra       -0.083     -0.097      0.124      0.827     -0.165      0.274
pmdec1          9.286     -6.037      6.882
pmdec2          5.316     -1.194      1.139
pmdec3          3.106     -0.572      0.568
pmdec4          2.512     -0.327      0.340
pmdec5          1.133     -1.909      1.306
lsigpmdec      -0.144     -0.111      0.132      0.717     -0.161      0.254
bv            209.129    -51.321     77.288
lsigbv          2.177     -0.031      0.040    150.286    -10.398     14.405
bfeh           -1.485     -0.025      0.029
lsigbfeh       -0.244     -0.014      0.015      0.570     -0.018      0.020
bpmra          16.592     -3.052      2.329
lsigbpmra       1.011     -0.030      0.021     10.261     -0.684      0.515
bpmdec         -5.830     -0.118      0.107
lsigbpmdec      0.650     -0.008      0.008      4.469     -0.081      0.087

Figure 17¶

InĀ [111]:
importlib.reload(st)

# Use MCMC results for plotting
meds = mcmc.meds # Store MCMC results in MCMeta
original_params = MCMeta.initial_params.copy()

# Convert MCMC results to the format expected by the plotting function
mcmc_for_plotting = {}

# Convert individual spline points back to arrays
vgsr_points = np.array([meds['vgsr1'], meds['vgsr2'], meds['vgsr3'], meds['vgsr4'], meds['vgsr5']])
pmra_points = np.array([meds['pmra1'], meds['pmra2'], meds['pmra3'], meds['pmra4'], meds['pmra5']])
pmdec_points = np.array([meds['pmdec1'], meds['pmdec2'], meds['pmdec3'], meds['pmdec4'], meds['pmdec5']])

mcmc_for_plotting = {
	'vgsr_spline_points': vgsr_points,
	'lsigvgsr': meds['lsigvgsr'], 
	'feh1': meds['feh1'],
	'lsigfeh': meds['lsigfeh'],
	'pmra_spline_points': pmra_points,
	'lsigpmra': meds['lsigpmra'],
	'pmdec_spline_points': pmdec_points,
	'lsigpmdec': meds['lsigpmdec'],
	'bv': meds['bv'],
	'lsigbv': meds['lsigbv'],
	'bfeh': meds['bfeh'],
	'lsigbfeh': meds['lsigbfeh'],
	'bpmra': meds['bpmra'],
	'lsigbpmra': meds['lsigbpmra'],
	'bpmdec': meds['bpmdec'],
	'lsigbpmdec': meds['lsigbpmdec']
}

# Update MCMeta with MCMC results
MCMeta.initial_params = mcmc_for_plotting

print("\nPlotting MCMC-based Gaussian mixture...")
plt_mcmc = st.StreamPlotter(MCMeta)
plt_mcmc.gaussian_mixture_plot(showStream=False, background=True, show_model=True, show_total=True)

# Restore original initial_params
MCMeta.initial_params = original_params

#plt.suptitle('MCMC-based Gaussian Mixture Model', fontsize=14, y=0.98)
plt.tight_layout()
plt.show()
Plotting MCMC-based Gaussian mixture...
No description has been provided for this image

Membership Probability¶

We now have our final gaussian mixture model parameters! We can go ahead and then calculate the probability that each star belongs to the 'stream' population.

Figure 18¶

InĀ [112]:
# Calculate membership probabilities using mcmc.memprob function
stream_prob = mcmc.memprob()


# Create membership probability histogram
fig, ax = plt.subplots(1, 1, figsize=(7, 6))

#-----------------------#
min_prob = 0.5  # Change this to adjust the stream probability threshold
#-----------------------#

ax.hist(stream_prob, bins=50, alpha=0.7, color='blue', 
        label=f'All Stars ({len(stream_prob)} total)', density=True)

ax.axvline(0.5, color='black', linestyle='--', alpha=0.7, label='50% threshold')

# Format plot
ax.set_xlabel('Membership Probability', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title(f'Membership Probability Distribution for {streamName}', fontsize=14)
ax.legend()
ax.grid(alpha=0.3)
stream_funcs.plot_form(ax)

# log yaxis
ax.set_yscale('log')

plt.tight_layout()
plt.show()

# Print statistics for different probability thresholds
thresholds = [0.3, 0.5, 0.9]
print("\nMembership statistics:")
print("Threshold | Count")
print("----------|-------")
for thresh in thresholds:
    count = len(stream_prob[stream_prob >= thresh])
    print(f"   ≄{thresh:3.1f}   | {count:5d} ")
Calculated membership probabilities for 1666 stars
Membership probabilities range from 0.000 to 1.000
Mean membership probability: 0.012
Stars with >50% probability: 21
Stars with >70% probability: 20
Stars with >90% probability: 19
No description has been provided for this image
Membership statistics:
Threshold | Count
----------|-------
   ≄0.3   |    21 
   ≄0.5   |    21 
   ≄0.9   |    19 

Now lets take a look at the stream you found!

InĀ [113]:
importlib.reload(st)

plt_enhanced = st.StreamPlotter(MCMeta) 
plt_enhanced.plot_params['sf_in_desi']['alpha'] = 0.5
plt_enhanced.plot_params['sf_in_desi']['zorder'] = 10
plt_enhanced.plot_params['background']['alpha'] = 0.2
plt_enhanced.plot_params['background']['s'] = 1
plt_enhanced.plot_params['optimized_spline_line']['lw'] = 1
plt_enhanced.plot_params['initial_spline_line']['lw'] = 1

fig, ax = plt_enhanced.sixD_plot(
    showStream=False, 
    show_sf_only=False, 
    background=True,
    show_initial_splines=True,
    show_optimized_splines=True,
    show_sf_errors=False,
    show_mcmc_splines=True,
    show_membership_prob=True, stream_prob=stream_prob, min_prob=min_prob, show_residuals=False
)

#plt.tight_layout()
plt.show()
No description has been provided for this image
InĀ [114]:
importlib.reload(st)

plt_enhanced = st.StreamPlotter(MCMeta) 
plt_enhanced.plot_params['sf_in_desi']['alpha'] = 0.5
plt_enhanced.plot_params['sf_in_desi']['zorder'] = 10
plt_enhanced.plot_params['background']['alpha'] = 0.2
plt_enhanced.plot_params['background']['s'] = 1
plt_enhanced.plot_params['optimized_spline_line']['lw'] = 1
plt_enhanced.plot_params['initial_spline_line']['lw'] = 1

fig, ax = plt_enhanced.sixD_plot(
    showStream=False, 
    show_sf_only=False, 
    background=False,
    show_initial_splines=False,
    show_optimized_splines=False,
    show_sf_errors=False,
    show_mcmc_splines=True,
    show_membership_prob=True, stream_prob=stream_prob, min_prob=min_prob, show_residuals=True
)

#plt.tight_layout()
plt.show()
No description has been provided for this image
InĀ [115]:
mcmc.save_run()
Calculated membership probabilities for 1666 stars
Membership probabilities range from 0.000 to 1.000
Mean membership probability: 0.012
Stars with >50% probability: 21
Stars with >70% probability: 20
Stars with >90% probability: 19
Saved 21 high-probability members to: ./runs/Gaia-9-I21_250811_test/Gaia-9-I21_phi2_spline_50%_mem.fits
Saved 1666 total stars to: ./runs/Gaia-9-I21_250811_test/Gaia-9-I21_phi2_spline_all%_mem.fits